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The totally asymmetric simple exclusion process in discrete time is considered on finite rings 
with fixed number of particles. A translation-invariant version of the backward-ordered sequential 
update is defined for periodic boundary conditions. We prove that the so defined update leads to 
a stationary state in which all possible particle configurations have equal probabilities. Using the 
exact analytical expression for the propagator, we find the generating function for the conditional 
probabilities, average velocity and diffusion constant at all stages of evolution. An exact and explicit 
expression for the stationary velocity of TASEP on rings of arbitrary size and particle filling is 
derived. The evolution of small systems towards a steady state is clearly demonstrated. Considering 
the generating function as a partition function of a thermodynamic system, we study its zeros in 
planes of complex fugacities. At long enough times, the patterns of zeroes for rings with increasing 
size provide evidence for a transition of the associated two-dimensional lattice paths model into a 
clustered phase at low fugacities. 
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I. INTRODUCTION 

The totally asymmetric simple exclusion process (TASEP) is, probably, the simplest model in the kinetic theory of 
many interacting particles Being one of the few exactly solvable models of non-equilibrium statistical mechanics, 

the TASEP describes a nontrivial evolution of the system from an initial configuration to a steadystate. While the 
long-time stage of this evolution and the steady state itself have been studied in details 0jI1jIEHj0IEI3) much less 
is known about the initial short-time stage even for relatively simple small systems. 

The exact solution of the master equation obtained by Schiitz [l(J gives a complete description of the process for 
arbitrary time intervals. However, the considered case of finite number of particles on an infinite chain has a trivial 
steady state in which all particles become noninteracting. By this reason, the Schiitz's solution was generalized for 
the ring geometry where a steady state of finite density exists. The ring solution is valid for both continuous and 
discrete time dynamics. In this paper we use the latter to demonstrate how the system tends to its stationary state. 

The discrete time TASEP, to be solvable on a ring, needs a special translation-invariant version of the backward- 
ordered sequential dynamics. The backward-ordered update is well known in both particle-oriented and site-oriented 
versions |l2lll3| . however, its original definitions break the invariance of the dynamics either with respect to transla- 
tions or cyclic permutations. That is why we begin with an invariant re-definition of the site-oriented backward-ordered 
sequential update for periodic boundary conditions. Then we show that our update leads to a stationary state where 
all possible particle configurations have equal probabilities. The rest of the paper is purely computational: using the 
exact formula for the conditional probabilities we find the corresponding generating function, which is an analog of 
the partition function in equilibrium statistical mechanics. From it we obtain the evolution of the average velocity 
and diffusion constant. The zeros of the generating function present a geometric portrait of the process and give 
evidence for a clustering phase transition in an associated two-dimensional lattice path problem. 
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II. THE MODEL 



Since details of the particle dynamics may have a strong effect on the time evolution of the system, we pay due 
attention to the definition of an invariant backward-ordered sequential update for the discrete-time TASEP on a ring. 
We consider a finite ring of L sites, each of which may be empty or occupied by a particle. In TASEP particles are 
allowed to hop only to the right: if site i is occupied and site i + 1 is empty, the particle may hop to site i + 1 with 
probability z, or remain on site i with probability 1 — z. In the case of a site-oriented discrete-time dynamics, one 
time-step corresponds to updating each site of the ring in a given order. According to the standard definition of the 
site-oriented backward-ordered sequential dynamics [lJl, the labeling of the sites is fixed, % = 1, 2, . . . , L, (L + i = i), 
and the local updates are applied in the order L, L — 1, . . . , 1. It is readily seen that the fixed beginning of this sequence 
breaks the translation invariance of the system: if site L is occupied by a particle, the possible updates depend on 
whether site 1 is empty or not. The same holds true for the particle-oriented backward-ordered sequential dynamics 
|12| : since particles cannot overtake each other, their order fx = 1, 2, . . . , M is preserved in time and the local updates 
always begin with particle M. Hence, the possible outcomes depend on the presence of empty sites between particles 
fj, = 1 and [i = M . In the case of a factorized steady state this feature leads to a single-particle weight for fi = M 
different from the (equal) weights of the other particles • 

Our invariant, site-oriented backward-ordered sequential update is defined as follows: at each time step one chooses 
an empty site, say i, and begins the sequence of local updates in the order i, i — 1, . . . , i + 1. Obviously, the choice of 
the initial empty site does not affect the possible outcomes of the updates in each time step. 

A distinctive property of the so defined update is the invariance with respect to the simultaneous reversion of space 
and time (this is not the case for the parallel update). To prove this, consider a time step t — > i + 1 from an arbitrary 
configuration C to configuration C allowed by the backward sequential dynamics. A particle jumps from site x to 
empty site x + 1 with probability z, and stays at the same site x with probability (1 — z) if x + 1 is empty, and with 
probability 1 otherwise. Reversing both the time and space directions, we obtain a step t + 1 — ► t from C to C . 
According to the reversed dynamics, the step from x + 1 to x has probability z\ immobility at x + 1 has probability 
(1 — z) if x is empty and 1 otherwise. Noticing that the number of factors (1 — z) is equal to the number of clusters 
of neighboring immobile particles, we see that the total transition probabilities for the direct and reversed dynamics 
coincide. 

Consider now the equiprobable distribution Pq(C) = (p) at the initial moment of time. The probability of a 
configuration C at the moment T, Pt(C ) is given by 

Pt(C')=( L ) ^2W T (C^C'), (1) 
^ ' c 

where Wt(C — > C ) is the transition probability. Continuing Eq. Q and using the reversed dynamics, we can write 

W) = (£) 1 J2 w t( c '^ c ) = (p) (2) 

c 

due to the normalization identity 

^W T (C' ->C) = 1. (3) 

c 

Therefore, if -Po(C) is a stationary equiprobable distribution, it remains equiprobable after T time steps. The properties 
of the transition probabilities guarantee its existence and uniqueness. Thus, we have proved that under the above 
defined backward-ordered stochastic dynamics, the system of fixed number of particles on a finite ring, starting from 
arbitrary initial configurations, evolves towards the equiprobable stationary state. 
A basic object for our consideration is the conditional probability 

Pt(C\C ) = Pt(x 1 ,...,Xp\x 1 ,...,x° p ) (4) 

to find P particles on lattice sites < x\ < x-i < ... < xp < L after T time steps from the initial state < x\ < x® < 
... < x P < L. The Bethe ansatz decouples a complicated many-particle dynamics into simple single-particle Bernoulli 
processes characterized by the probability 



B(N,T)=[ „ )z«(l-zy-» (5) 
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of N advances for T time steps. Here, we assume that B(N,T) = for N > T. The solution obtained in ^lj reads 

oo oo 

P T (C\C a )= - J2 (-l) (p - 1)E -i"*detM, (6) 

tlx = — OO Tip — — OO 

where the elements of the P x P matrix M are 

M ij = F Sij (x° i ,x j +n j L\T) 

with 

s ij = ( p - l ) n j - ^2 n k + j - i, 



and the functions F m (x®, Xj\T) are defined as 

F m {xl Xj \T) = £ ( fc + m 7 1 ) B(x° - k, Xj \T), 



fc=0 

if integer m > 0, and 



F m {xl xj\T) = fj(-l) fc C ™ ) - k, Xj \T), 
fc=0 ^ ' 



if integer m < 0. It should be noted that the infinite sums in © are actually bounded from above and below for finite 
T, since the arguments of non-zero Bernoulli function are restricted by the condition N < T. 



III. GENERATING FUNCTION 

The probability Pt(C\C°) involves all transitions from C° to C with arbitrary numbers of rotations ni, ri2, ■ ■ ■ , rip, 
of the particles around the ring for time T. It is convenient to select explicitly the transitions in which the total 
number of steps Y of all particles is fixed: 

Y = L(ni + ... + np) + (xi + ... + x P )-(x1 + ... + x%). (7) 
The corresponding probability has the form : 

P%(C\C°)= Yl (-l) {P ' 1)j: ^ ini detM, (8) 

{m,—,np}' 

where the summation over the set {ni, . . . , np} is restricted by condition J7J. 
The generating function for Prf (C\C°), 

oo 

P^{C\C°) = w Y P£{C\C°), (9) 

Y=0 

can be written as a unrestricted weighted sum over the numbers of rotations 

oo oo 

P^(C\C")= w y (-l) (p - 1)S -i"'detM (10) 

n\ — — oo np — — oo 

Summation over all possible configurations C gives a generating function At(w,z) which can be considered as a 
"canonical partition function" of the discrete process for time T: 

K t {w,z) = Y,Pt{C\C\ (11) 

c 
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Actually, such a generating function has been introduced by Derrida and Lebowitz 9] in a study of the large deviation 
properties of the continuous-time TASEP on a ring. 

An example of At(w, z) for the case L = 5, P = 3 , time T = 3 and initial positions x\ = 0, x® = 1, x% = 2 is 



A 3 (w,z) = l-3z + 3wz + 3z 2 - I2wz 2 + 9w 2 z 2 - z 3 + 19wz 3 -36w 2 z 3 + 18w 3 z 3 - I5wz 4 
+ 57 w 2 z 4 - 63 w 3 z 4 + 21 w 4 z 4 + 6 w z 5 - 45 w 2 z 5 + 87 w 3 z 5 - 69 w 4 z 5 + 21 w 5 z 5 - w z 6 
+ 18 w 2 z 6 - 60 w 3 z e + 85 w 4 z 6 - 60 w 5 z 6 + 18 w 6 z 6 - 3 w 2 z 7 + 21 w 3 z 7 - 48 w 4 z 7 



60 w 5 z 7 - 39 w 6 z 7 + 9 w 7 z 7 - 3 w 3 z 8 + 12 w 4 z 8 - 24 w 5 z 8 + 27 w 6 z 8 - 15 w 7 z 8 



,.9 ~9 



+ 3 w* z 8 - w 4 z y + 3 w° z y - 6 w b z J + 6 w ' z y - 3 w* z 9 + w M z 

Larger polynomials can be obtained from Eqs. l(6"|l- i|ll|) by using the package MATHEMATICA. We point out that 
the degree of the polynomials considered here reaches 160, and the number of digids in their integer coefficients is up 
to 70. 

Note that due to the probability normalization condition, At(w = 1, z) = 1 for all values of the model parameters. 
The negative terms in the expression for At(w, z) appear due to the factors (1 — z) attached to each cluster of immobile 
particles at each time step. Using (fill) , one easily obtains the average distance travelled by all particles for time T: 



00 = 



dA T (w) 



dw 



the average velocity 



and diffusion constant 



P T ' 



A = 



1 <9 2 A T («;) 



T dw 2 



(Y) - {Yf 
T 



(12) 



(13) 



(14) 



The stationary velocity can be obtained from At(w, z) for the one-step evolution if one takes all the initial configu- 
rations with equal weight (p) : 



c° c 



p \p / 

The above expression can be cast, see the Appendix, in the following explicit form: 



(15) 



v st {L,P,z) 



(P-1)\(L-P) A 



(L-iy. 



E 

i=l 



L-P 
L-l 



2 F 1 (1,1-P,2-L;z), 



(16) 



where 2F1 is the hypergeometric function. Particularly, in the thermodynamical limit (infinite chain with fixed density 
of particles p — P/L) we recover the known result [l3| 



pz 



On the other hand, the continuous time limit yields the result obtained by Derrida and Lebowitz 

L-P 



L - 1 



(17) 



(18) 



Several values of v s t for small lattices (L, P) at z = 0.5 are 



«*(4,2) = - 



^(5,2) = ^; M5,3) = g 
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w s t(6,2) = — ; v st (6,3) = — ; v st (6,4) = — 

Convergence of velocities to their stationary values can be obtained from (|13fl . Fig. ^ shows how the exactly 
calculated velocities for the case L = 8, P — 4, z — 0.5 and time intervals T = 10, 20, 30, 40 tend to the stationary 
value v st (8,4) = 209/560. 

We can illustrate also the continuous-time behavior on large time scales by rescaling the variables z — > 0, T — > oo, 
Tz -» t, see Fig. 
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FIG. 1: Convergence of the finite-time velocity to the stationary value for L = 8, P — 4. 





FIG. 2: The tendency of current j = Pv (a) and diffusion constant (b) to stationary values (the horizontal line) for L = 5, 
P — 3 in continuous time. The two branches correspond to different initial conditions. 

As it is clearly seen, the velocity and the diffusion coefficient relax to their known stationary values, thus giving 
evidence of the correctness of the general expressions derived here. Although the concrete way of approach to 
stationarity depends strongly on the initial condition, the long-time relaxation seems to be exponential, at least for 
small systems. 



IV. ZEROS OF THE PARTITION FUNCTION 



^,From the equilibrium theory by Yang and Lee |14j it is known that the zeros of the partition function of a finite-size 
system in the complex fugacity plane approach the real axis when the system size increases. By calculating the line of 
zeros in the thermodynamic limit, and their density near the real axis, one can exactly locate the transition point and 
obtain the order of the phase transition. Recently, essential progress has been made in the extension of the Lee- Yang 
theory to the non-equilibrium transitions between stationary states of the TASEP. Blythe and Evans , 0] , found 
that the distribution of the zeros of the normalization factor for the continuous-time TASEP with open boundaries in 
the plane of complex injection rate agrees with the predictions of the Lee- Yang theory, see also |17| . Exact mappings 
of the normalization factor for the discrete-time parallel update onto several equivalent two-dimensional lattice path 
problems were obtained in [Djj. The applicability of the concepts of finite-size scaling and universality with respect 
to the different types of updates for the open TASEP has been established in 19]. For an overview on the recent 
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breakthroughs in the extension of the Lee- Yang theory to non-equilibrium phase transitions we refer the reader to 
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FIG. 3: Two configurations of paths are shown. The ovals surround different clusters of particles standing during tree time- 
steps: in part (a) the cluster of four standing particles contributes a total weight of y 3 , and in part (b) the two clusters, of two 
standing particles each, contribute a total weight of xy 6 . 

Let us now re-interpret the generating function (jl 1|> as a partition function of a two-dimensional lattice path 
problem, in which the horizontal coordinates correspond to the original lattice sites, and the vertical direction is the 
discrete time. Each step to the right has a statistical weight x — wz, and each vertical step has a weight y = (1 — z) 
if the target site for the next step is empty and the weight 1 if it is occupied. Fig. [3] illustrates the equivalent 
two-dimensional path problem. The statistical weights x, y, 1 are ascribed to all steps of the paths in accordance with 
the TASEP dynamics. In terms of these new variables, takes the form of a polynomial Zx{x,y) with positive 
coefficients. For instance, the polynomial corresponding to A^lw^z) is 

Z 3 (x, y) = x 9 + 3x 8 y + 3x 7 y{l + 2y) + 3x 6 y{l + 3y + 2y 2 ) + 3x 5 y 2 (2 + 4y + y 2 ) + x 4 y 2 (3 + Wy + 7y 2 + y 3 ) 

+ 3x 3 y 2 {l + 2y + 2y 2 + y 3 ) + 3x 2 y 3 (l + y + y 2 )+ xy 3 (l + y + y 2 ) + y 3 . (19) 

Obviously, the probability normalization condition now implies Zt(x,1 — x) = 1. Inspired by the success of the 
Lee- Yang theory in explaining how singularities may build up in the thermodynamic limit, we attempt to clarify the 
analytical structure of Zt{x, y) by studying the location of its zeros in the planes of complex fugacities x and y, for 
sufficiently long times T, different ring sizes L at fixed particle density P/L and given initial conditions. The patterns 
of zeros obtained in the complex- a; and complex-y planes are illustrated in Figs. 0] and 03 respectively, for three cases: 
L = 6, P = 3, T = 40; L = 8, P = 4, T = 40; and L = 10, P = 5, T = 30. 
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FIG. 4: The pattern of zeros of the partition function Zt(x, 0.5) in the complex- a; plane for L — 6, P — 3, T = 40 (solid 
squares), L — 8, P = 4, T — 40 (open circles), and L — 10, P — 5, T = 30 (open triangles). Note that the number of branches 
equals the number of particles P. 
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FIG. 5: The pattern of zeros of the partition function Zt(0.5,j/) in the complex-?/ plane for L = 6, P — 3, T = 40 (solid 
squares), L = 8, P = 4, T = 40 (open circles), and L — 10, P — 5, T = 30 (open triangles). 



By inspecting these figures, one may speculate on the tendency of the zeros to approach the positive fugacity axis 
with increasing L at fixed P/L. That may signal the appearance of a kind of condensation phase transition in the 
related thermodynamic model with fugacities x and y. Arguments in favor of such a conjecture are given in the next 
section. 



V. CONDENSATION 



To reveal the nature of the possible phase transition, it is convenient to use the mapping of the TASEP onto the 
zero-range process (ZRP), see [2lJ for a recent review. Consider a ZRP each site of which is identified, in the same 
sequential order, with one of the L — P empty sites of the TASEP. Let the number of particles at each site of the ZRP 
be equal to the number of neighboring occupied sites of the TASEP just on the left of the corresponding empty site. 

The interpretation of the partition function Zt{x, y) in terms of the ZRP is the following. Each jump of a particle 
from site i to site i + 1 [i = 1, 2, . . . , L — P; i + L — P = i) for the unit time step from t to t + 1 carries a weight x, 
as in the TASEP. Each site occupied by least one particle which does not move during that time step has the weight 
y. The total contribution to ZT(x,y) from a particular realization of the process is the product of weights over all 
time steps from to T. Comparing the contributions with different factors y n and y m , n < m, and equal number of 
spatial steps (i.e. factors x k ), we see that the decrease in y leads to reduction in the average number of sites occupied 
by immobile particles and, therefore, to an effective clustering attraction between them. For sufficiently small y, i.e. 
sufficiently strong attraction, a condensation of particles may take place at a site (or at a number of sites). 

Real-space condensation phenomena are known to take place in the stationary states of both the heterogeneous 
ZRP, with site-dependent hopping probabilit ies, and the homogeneous ZRP, where the hopping rates depend only on 
the occupation number of the departure site [2l| . In the former case the P particles of the ZRP are viewed as bosons 
which may have N = L — P energies determined by the different single-site hopping rates; the energy levels are defined 
so that the ground state corresponds to the lowest hopping rate. In the discrete-time TASEP with quenched random 
hopping probabilities of the particles |12|. the phase transition takes place between a low-density inhomogeneous 
phase, with a traffic jam behind the slowest particle (and empty space in front of it), and a high-density congested 
phase. 

In our case the condensation is due to a specific interplay between the immobility probability y and the fugacity 
x associated with the hopping probability. Turning back to the effective two-dimensional path problem, we see that 
vertical bonds, corresponding to immobile neighboring particles, tend to form compact clusters, or jams in the TASEP 
traffic. This situation is shown in the ovals in Fig. [3] The cluster of four standing particles in Fig. |3| (a) lives three 
time steps and brings the total weight y 3 . The configuration of paths in Fig. |3| (b) corresponds to two clusters 
containing two standing particles each. The weight of the configuration in the second oval is xy 6 that is much less 
than y 3 for small y. Varying parameters x, y, one may attain an arbitrary strong effective attraction between clusters. 



VI. DISCUSSION 



Here we considered the exact evolution of the TASEP on small rings during finite intervals of discrete time. The 
results for the velocity and diffusion constant are not surprising. As expected, both of them tend to their stationary 
values exponentially, although the particular features of the convergence depend on the initial conditions. 



8 



A more interesting object is the generating function of the process. Considering it as a function of two independent 
variables, fugacities of immobile particles and particle jumps, we studied the loci of its zeros in the planes of complex 
fugacities. When the location of the zeros in the right-hand half-plane is well behaved, a finite-size scaling analysis can 
provide important results concerning the very existence, the location, and the characteristics of the conjectured phase 
transition. However, the study of larger lattices and longer evolution times poses serious computational problems and 
will be presented in a forthcoming work. Although the tractable system sizes are too small, the obtained patterns 
suggest that with increasing the lattice size at a fixed particle density, the zeros may concentrate at a complex phase 
boundary which cuts the positive real axis at the transition point. When the corresponding fugacity crosses that 
point, a phase transition in the related two-dimensional equilibrium lattice path model should occur. We interpret 
this transition as a condensation of immobile particles. However, several problems remain unsolved. 

The first one is to prove the existence of a phase transition in the framework of equilibrium statistical mechanics. To 
this end, one has to define an order parameter and a low-temperature phase which is destroyed when the temperature 
increases. If one succeeds in defining an interface between the phases, the proof can be based on contour arguments 
which compare the entropy and energy of the boundary of the low-temperature clusters. 

The second problem is the kind of the phase transition. According to the Lee- Yang theory, it can be determined 
by evaluating the asymptotic density of zeros around the real axis. However, the fulfillment of this program needs 
consideration of much larger systems than the ones considered here. 

It is of interest also to construct the complete phase diagram of the model in the (z,w) plane, or in the (x,y) 
plane, and to clarify the relationship between the dynamics of the TASEP and the thermodynamics of the effective 
two-dimensional model. 
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Appendix 

Here we present the derivation of expression (|16f) for the stationary velocity in finite-size systems. 
After summation over all final configurations C, Eq. I|15|l can be written as a sum over all clusters (connected 
sequences of neighboring particles) in the initial configuration C°: 

v ' C° i=l 

where m(C°) is the number of clusters in configuration C° and Ni(C°) is the number of particles in the i-th cluster 
of the configuration C . For further calculations, it is convenient to introduce the following generating function of 
the stationary velocity: 

T(L,z,s)=Y, P ( p )vst(L,P,z) S p . (21) 
p=i ^ ' 

It is readily seen that the stationary TASEP can be mapped on the Fortuin-Kasteleyn representation of the one- 
dimensional Potts model with external field [24[ • To this end, we introduce the ring graph G with L vertices located 
at the centers of the bonds of the TASEP ring R. Let G' be a subgraph of G, such that two neighboring vertices are 
connected by an edge if and only if the common site of the corresponding bonds of R is occupied. Thus, we obtain a 
one-to-one correspondence of the TASEP configurations on R and the subgraphs of G. Note that the graph G itself 
and the edgeless subgraph correspond to the fully occupied (P = L) and the empty (P — 0) rings, respectively. The 
velocity generating function (|21|l can be represented as a sum over all subgraphs: 

m z,s) = yzt z E * e(G ' } E (1 - ^ (GVl ) - (1 - ^) A (22) 
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where e(G') and k(G') are correspondingly the number of edges and the number of connected components of G', 
rii{G') is the number of vertices in the i-th connected component of G' . The last therm is a correction to the fully 
occupied case. The partition function of the g-state Potts model in external field can be expanded over the spanning 
subgraphs as follows: 



KG') 



z Potts ( q ,v,w)= j2 ve(G,) n [(^+ir i(G,) +9-i 



G'CG i=l 

Here v and w are the coupling and external field parameters, respectively. If we denote 

1 d Zp tts 



f{v,w) 



(w + l) L dq 



9=1 



then the generating function l|22[l can be written as 



r(L, *,*) = --- f(s,0)--^—f(s,—)- ' 



1 — z 1 — z z 1 — z 

By using the exact expression for the Potts model partition function, 

Z Potts (q, v, w) = {q- 2)v L + (A+) L + (A~) L , 

with 

l r 



{l-z L - l )s L . 



X ± = 



1 + 2v + w + vw ± y/ (q + 2v + w + vw) 2 — 4v(q + v)(l + w) 



after some simplifications we find 



T(L, z, s) = Lzs 



(1 + s)^" 1 - (zs) 



L-l 



1 + s — zs 

With the aid of the above exact expression we obtain the explicit expression for the stationary velocity. 
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